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ABSTRACT 

Time-dependent  reliability  is  the  probability  that  a 
system  will  perform  its  intended  function  successfully 
for  a  specified  time.  Unless  many  and  often  unrealistic 
assumptions  are  made,  the  accuracy  and  efficiency  of 
time-dependent  reliability  estimation  are  major  issues 
which  may  limit  its  practicality.  Monte  Carlo 
simulation  (MCS)  is  accurate  and  easy  to  use  but  it  is 
computationally  prohibitive  for  high  dimensional,  long 
duration,  time-dependent  (dynamic)  systems  with  a  low 
failure  probability.  This  work  addresses  systems  with 
random  parameters  excited  by  stochastic  processes. 
Their  response  is  calculated  by  time  integrating  a  set  of 
differential  equations  at  discrete  times.  The  limit  state 
functions  are  therefore,  explicit  in  time  and  depend  on 
time-invariant  random  variables  and  time -dependent 
stochastic  processes.  We  present  an  improved  subset 
simulation  with  splitting  approach  by  partitioning  the 
original  high  dimensional  random  process  into  a  series 
of  correlated,  short  duration,  low  dimensional  random 
processes.  Subset  simulation  reduces  the  computational 
cost  by  introducing  appropriate  intermediate  failure 
sub-domains  to  express  the  low  failure  probability  as  a 
product  of  larger  conditional  failure  probabilities. 
Splitting  is  an  efficient  sampling  method  to  estimate  the 
conditional  probabilities.  The  proposed  subset 
simulation  with  splitting  not  only  estimates  the  time- 
dependent  probability  of  failure  at  a  given  time  but  also 
estimates  the  cumulative  distribution  function  up  to  that 
time  with  approximately  the  same  cost.  A  vibration 
example  involving  a  vehicle  on  a  stochastic  road 
demonstrates  the  advantages  of  the  proposed  approach. 

1.  INTRODUCTION 

Reliability  is  an  important  engineering  requirement 
for  consistently  delivering  acceptable  product 
performance  through  time.  As  time  progresses,  the 


product  may  fail  due  to  time -dependent  operating 
conditions  and  material  properties,  component 
degradation,  etc.  The  reliability  degradation  with  time 
may  increase  the  lifecycle  cost  due  to  potential 
warranty  costs,  repairs  and  loss  of  market  share. 

Reliability  is  the  probability  that  the  system  will 
perform  its  intended  function  successfully  for  a 
specified  interval  of  time,  under  stated  operating  and 
environmental  conditions.  It  is  therefore,  related  to 
product  functionality  over  time  which  is  determined  by 
the  so-called  “hard”  and  “soft”  failures  [1],  In  a  hard 
failure  the  system  loses  functionality  due  to  a  complete 
breakdown  of  one  or  more  of  its  components,  while  in  a 
soft  failure  the  system  is  functional  but  one  or  more 
performance  measures  are  out  of  conformance.  The 
reliability  associated  with  the  hard  failure  is  important 
for  non-repairable  systems  where  the  replacement  or 
repair  of  a  failed  component  is  not  possible  and  the 
failed  system  is  removed  from  the  population.  In 
contrast,  repairable  systems  [2]  consist  of  multiple 
components  which  can  be  repaired  or  replaced  if  failed 
keeping  therefore,  the  original  system  in  the  population. 
In  this  research,  we  use  time-dependent  reliability 
concepts  associated  with  the  so-called  first-passage  of 
non-repairable  systems.  Among  its  many  applications, 
the  approach  can  be  used  to  reduce  the  lifecycle  cost  [3, 
4]  or  to  set  a  schedule  for  preventive  condition-based 
maintenance  [5]. 

The  time -dependent  probability  of  failure  (see  Eq. 
5  for  definition),  also  known  as  cumulative  probability 
of  failure  [3,  6],  is  calculated  by  the  following  exact 
relation  using  the  failure  rate 


Pf(0,7’)  =  l-(l-Pfi(0))exp|-{4)dfJ ,  (1) 
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where  Pf'(0)  is  the  instantaneous  probability  of  failure 
at  the  initial  time.  In  the  commonly  used  out-crossing 
rate  approach,  the  failure  rate 

Pit  <Tf  <  t  +  dr  I  T{  >t) 


Mt)  = 


lim- 

df— >0 


dr 


(2) 


where  Tf  is  the  time  to  failure,  is  approximated  by  the 

up-crossing  rate  [7] 

v+(r) 

_  p[g(X(t),t)-Si  <0og(X(r  +  Ar),r  +  Ar)-.S,  >o]  , 

Ar-»0,  At 

Ar>0 

(3) 

(see  Equations  5  and  6  for  notation)  under  the 
assumptions  that  1)  the  probability  of  having  two  or 
more  out-crossings  in  [r,r  +  Ar]  is  negligible  compared 
to  the  probability  of  having  exactly  one  out-crossing,  2) 
At  is  sufficiently  small,  and  3)  the  out-crossings  in 
[r,r  +  Ar]  are  statistically  independent  of  the  previous 
out-crossings  in  [0,r],  In  this  case,  the  number  of  up- 
crossings  N+(t,a)  for  a  threshold  a  is  a  Poisson 
process. 

The  out-crossing  rate  approach  was  first  introduced 
by  Rice  [8]  followed  by  extensive  studies  [6,  9-13] 
under  the  assumption  of  the  out-crossings  being 
statistically  independent  and  Poisson  distributed. 
Hagen  and  Tvedt  [7]  suggested  a  parallel  system 
reliability  formulation  to  compute  the  out-crossing  rate. 
It  uses  two  successive  time-invariant  analyses  based  on 
FORM,  and  the  binomial  cumulative  distribution  to 
calculate  the  probability  of  the  joint  event  in  Eq.  (3). 
This  approach  was  later  adopted  in  the  PHI2  method 
[6].  Methods  based  on  Poisson’s  distribution  and  the 
PHI2  method  compute  an  upper  bound  of  the 
probability  of  failure  of  Eq.  (3)  [14].  A  Monte-Carlo 
based  set  theory  approach  has  been  also  proposed  [1, 
15]  using  a  similar  approach  with  the  PHI2  method. 
Analytical  studies  such  as  in  [16,  17,  18]  revealed  that 
the  PHI2  based  approach  lacks  sufficient  accuracy  for 
non-monotonic  problems  such  as  vibratory  systems.  For 
this  reason,  analytical  approaches  were  proposed  in  [16, 
19,  20]  to  accurately  estimate  the  time -dependent 
probability  of  failure  considering  non-monotonic 
behavior. 

Although  the  out -crossing  rate  approach  can  easily 
estimate  the  time-dependent  probability  of  failure,  it  has 
two  potential  limitations.  First,  its  accuracy  may  be 
poor  because  of  the  Poisson  process  assumption  of 
independent  out-crossings  and  second,  it  may  require  a 
large  computational  effort.  An  analytical  FORM -based 
estimation  of  the  up-crossing  rate  (Eq.  3),  with  its  own 
accuracy  limitations,  must  be  performed  at  the  time 
instants  required  by  the  numerical  evaluation  of  the 
integral  in  Eq.  (1)  (e.g.  Gauss -Legendre  integration 
points).  If  the  probability  of  failure  /]((). 7’)  is 


calculated  at  different  times  T,  the  computational  effort 
increases  because  the  integration  points  change.  This 
can  increase  the  computational  effort.  The  first 
limitation  has  been  recently  improved  in  [18]  by 
considering  the  correlations  between  the  limit-state 
function  at  two  time  instants.  The  method  estimates  the 
up-crossing  rate  v'  it  )  by  solving  an  integral  equation 
involving  v+(t )  and  v++(t,tl ),  the  joint  probability  of 
up-crossings  in  times  t  and  tt  [21], 

This  paper  presents  a  simulation-based  method  to 
estimate  the  time-dependent  probability  of  failure  at 
different  time  instances.  Monte  Carlo  simulation  (MCS) 
can  handle  high-dimensional  problems,  and  general 
failure  definitions  allowing  us  to  handle  system 
reliability  problems.  However,  it  cannot  estimate 
efficiently  small  probabilities  because  the  number  of 
samples,  and  hence  the  number  of  system  analyses 
required  to  achieve  a  given  accuracy,  is  inversely 
proportional  to  the  failure  probability. 

Importance  sampling  techniques  [22]  are 
commonly  used  to  shift  the  underlying  distribution 
towards  the  failure  region  in  order  to  sample  rare  events 
more  efficiently.  They  require  however,  a  careful  choice 
of  the  importance  sampling  density  (ISD),  which 
requires  knowledge  of  the  failure  region.  For  low¬ 
dimensional  uncertain  systems  with  relative  simple 
failure  regions,  many  important  sampling  methods  have 
been  developed  (e.g.  [22,  23]).  However,  the 

application  of  importance  sampling  to  general  dynamic 
reliability  problems  where  the  random  excitation  is 
represented  by  a  large  number  of  discrete  random 
variables  is  still  an  active  research  area  with  limited 
practicality  for  such  problems. 

A  MCS  approach  was  proposed  in  [24]  to  estimate 
the  time-dependent  failure  rate  over  the  product 
lifecycle.  The  efficiency  of  the  method  was  further 
improved  using  an  importance  sampling  method  with  a 
decorrelation  length  [25]  in  order  to  reduce  the  high 
dimensionality  of  the  problem. 

Subset  simulation  [26,  27]  has  been  recently 
developed  as  an  efficient  simulation  method  for 
computing  small  failure  probabilities  for  general 
reliability  problems.  Its  efficiency  comes  from 
introducing  appropriate  intermediate  failure  sub- 
domains  to  express  the  low  probability  of  failure  as  a 
product  of  larger  conditional  failure  probabilities  which 
are  estimated  with  much  less  computational  effort.  In 
doing  so,  the  probability  of  a  rare  event  in  the  original 
probability  space,  is  replaced  by  a  sequence  of 
probabilities  of  more  frequent  events  in  conditional 
probability  spaces.  Because  it  is  very  challenging  to 
generate  samples  in  the  conditional  spaces,  subset 
simulation  with  Markov  Chain  Monte  Carlo 
(SS/MCMC)  [28,  29]  and  subset  simulation  with 
splitting  (SS/S)  [30-32]  methods  have  been  introduced. 


2 


In  this  paper,  we  propose  an  improved  subset 
simulation  method  with  splitting  by  partitioning  the 
original  high  dimensional  random  process  into  a  series 
of  correlated,  short  duration,  low  dimensional  random 
processes.  The  proposed  method,  called  subset 
simulation  with  splitting  and  partitioning  in  time 
(SS/SPT),  not  only  estimates  the  time -dependent 
probability  of  failure  at  a  given  time  (as  the  original 
SS/S  method  does)  but  also  estimates  the  failure  rate 
and  the  cumulative  distribution  function  (CDF)  up  to 
that  time  with  approximately  the  same  cost.  It  can  also 
handle  dynamic  systems  with  random  parameters  (e.g. 
random  stiffness  of  a  vehicle  suspension  system)  which 
the  original  SS/S  approach  does  not  consider. 
Estimation  of  the  CDF  is  very  important  for  many 
practical  problems  such  as  preventive  maintenance 
using  reliability  principles,  design  for  lifecycle  cost,  etc. 
Also,  metrics  such  as  the  failure  rate,  remaining  life, 
and  mean-time-to-failure  (MTTF)  which  are  functions 
of  the  CDF  of  time  to  failure  are  significant  for 
describing  the  system  reliability  and  performance 
characteristics. 

The  paper  is  organized  as  follows.  Section  2 
defines  the  dynamic  systems  we  consider  and  provides 
a  brief  introduction  of  the  problem  and  notation. 
Section  3  introduces  the  subset  simulation  approach  and 
the  existing  MCMC  and  splitting  sampling  schemes. 
Our  proposed  SS/STP  method  is  presented  in  Section  4, 
including  accuracy  bounds  and  computational  effort 
estimation.  Section  5  uses  an  example  of  a  quarter 
vehicle  with  uncertain  parameters  on  stochastic  terrain 
to  demonstrate  the  characteristics  and  advantages  of  the 
proposed  method.  Finally,  Section  6  summarizes, 
concludes  and  presents  future  work. 

2.  TIME  DEPENDENT  RELIABILITY 
EVALUATION 

We  consider  the  time-dependent  reliability  of 
dynamic  (rigid-body  or  vibratory)  systems  whose 
equations  of  motion  are  usually  discretized  in  time  and 
presented  in  a  state-space  form.  The  discretized 
equations  are  time  integrated  using  for  example,  a 
Runge-Kutta  method  or  Newmark-beta  method  [33], 
They  are  expressed  as 

X(t  +  At)  =  f(x(t),  Y,  U(r),  t) ,  (4) 

where  X(t  )g  41 ''  is  the  vector  of  uncertain  states 
xs  (t) ,  s  =  1,2,. . p  at  time  t.  At  is  the  integration  time 
step,  Y  g  9U  is  the  time-independent  vector  of  random 
variables  (e.g.  system  parameters  Ys  e  S.R''  and 
excitation  parameters  Ye  g  s.R''  ),  and 

U(f)=U(Y,f)G9T  is  the  time-dependent  vector  of 
excitation  random  processes  (e.g.  road  elevation  at  a 
vehicle  tire  location  through  time).  Both  X(r)  and  U(r) 


are  implicit  functions  of  Y .  The  trajectories 
X  =  {x(f),  t  g  [o,  /  ]]  of  all  states  (sample  functions  of 
corresponding  random  processes)  are  calculated  at 
discrete  time  instances  t  =  t ■,  j  =  0,l,...//stcp  where 

A/,L.p  is  the  number  of  time  integration  steps  over  the 
period  [0,  T\ .  For  illustration,  if  p  —  1 ,  q  =  2  and 
r  =  1,  Eq.  (4)  becomes 

x{(t  +  At))=  f(xl(t),yl,y2,u(t),t)  where  the  single  state 
x{  at  time  t  +  At  is  a  function  of  the  state  at  time  t,  y, 
and  y0  are  the  realizations  of  the  two  random  variables 
lj  and  Y2,  and  u(t)  is  the  value  of  a  sample  function  of 
the  random  process  u(t)  at  time  t.  Time  integration  of 
Eq.  (4)  provides  the  system  response 

*i  =  fa  (t),  t  g  [o,r]}. 

The  computational  effort  to  solve  for  the  system 
states  x(r  +  Ar)  at  time  t  +  At  as  a  function  of  the 
states  X(r)  at  time  t,  is  considered  one  function 
evaluation.  It  is  important  to  note  that  the  system  states 
X(t  +  Ar)  are  a  function  of  the  states  X(r) which  are  in 
turn  a  function  of  the  states  X(r-  Ar).  Thus,  in  order  to 
calculate  X(r),  we  need  the  solution  at  all  previous 
times.  For  this  reason,  counting  the  number  of  sample 
function  (or  trajectory)  evaluations  instead  of  the 
number  of  function  evaluations  is  preferable  for 
dynamic  systems.  This  is  different  from  problems  with 
explicit  in  time  limit  states  where  a  simple  function 
evaluation  can  be  performed  at  any  time  t. 

Our  goal  is  to  calculate  the  time-dependent 
probability  of  failure 

Pf  (o,  r)  =  p{3t  g  [o,  r] :  g{x(t),t)  >  St  } ,  (5) 

where  g:i Rp  —>9?  is  a  function  that  maps  X(r)  to  a 
response  of  interest  and  St  is  a  given  threshold  value. 
Because  of  the  time-dependent  uncertain  states  in  X(r), 
G(t)=g(x(t),t)-S,  is  a  random  process  which  can  be 
viewed  as  a  collection  of  random  variables  at  different 
time  instances  t.  Since  we  consider  a  first  excursion 
failure  problem  in  Eq.  (5),  the  failure  domain,  as  well  as 
an  event  therein,  can  be  defined  as 

F={  •  (6) 

The  system  operates  properly  and  is  called  safe 
if  g(x(f),t)<  St,Vt  g  [o,r].  The  system  is  considered 
failed  if  3t  g  [o,r]:  g(x(f),f)>  St  .  The  equation 
g(x(f),f)-.S’,  =0  is  the  time-dependent  limit  state 
surface  (LSS). 

We  consider  a  non-repairable  system  where  if 
g(x(?7),rj ,  )>  A,  for  r;  g  [o,  /  ] ,  the  system  fails  and  is 

removed  from  the  population  for  t>tJ.  In  this  case, 
Eq.  (5)  can  be  approximately  rewritten  as 
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Mo.r)=p|u*(xMf>st,  tjs[oA 

U=1  J  ,  (7) 

=  /’i  max  ,.-(\(r  l/>  St,  tj  e  [0,7']j 

where  r  is  the  number  of  load  random  variables 
generated  from  the  random  process 
G(t)  =  x(x(r\t)- S,  at  discrete  times  tf  ,  j  =  l,2,...,r  . 
For  convenience,  if  Na  is  the  number  of  time 
integration  steps,  we  let  r  =  /V9cp .  If  the  time  interval 
[o,r]  is  discretized  with  a  uniform  time  step  At ,  we 
have  7Vstep  =  .  According  to  Eq.  (7),  the  time- 

dependent  reliability  evaluation  can  be  transformed  into 
a  time -independent  system  reliability  evaluation 
involving  r  =  IVaep  correlated  random  variables.  If  Af 
is  small,  the  time-independent  system  reliability 
problem  is  of  very  high  dimension  since  Naep  is 
inversely  proportional  to  Af . 


3.  SUBSET  SIMULATION 


This  section  provides  a  brief  introduction  to  subset 
simulation  and  the  existing  MCMC  and  splitting 
sampling  schemes.  Given  the  failure  event  F  of  Eq.  (6), 
a  nested  sequence  of  m  failure  domains  is  formed  using 
a  set  of  increasing  threshold  levels 

.S',1  <  .S,1  <  ■ 

•  •  <  S'”  =  St .  The  failure  domains 

F‘  for 

i  =  1,2,..., 

m  are  defined  as 

F'=j  maxg(x(r),r)>  At'|  , 

[fE[o,r]  J 

(8) 

where 

^7!  — j  ^72  _ pm  _  p 

(9) 

and 

l 

Fl  =f)F‘ ,  I  =l,2,...,m . 

GO) 

i=i 


Figure  1  shows  pictorially  the  m  failure  domains 
F' ,  i  =  l,2,...,mfor  a  hypothetical  case  with  two 
random  variables. 

Using  the  definition  of  conditional  probability,  we 
have  [31] 


Figure  1.  Pictorial  representation  of  failure 
subdomains 

Pf  =  p(f)=  p(Fm)=  P 

v>=i 

f  m- 1  S  f  m- 1  A 

=  P  F'"  P)F;  P  f)F' 

v  i=l  J  V  i=l  J 


m— 1  /  \  (  m— 1 

=  n^r+iH^ 

i=m— 2  v  i=l 

,  (ID 

in — I  /  \  2  ^ 

i= 2  V  i=l  J 

=  p(f1)-]Qp(f'+1|f'') 

i=1 

m 

=  /Ml/r 

i-2 

m 

i= 1 

where  the  conditional  probability  , 

i  =  2,...,m  is  denoted  by  /]'  and  is  denoted  by 

Pf  .  Eq.  (11)  indicates  that  the  probability  of  failure  /] 
can  be  expressed  as  the  product  of  a  sequence  of  larger 
conditional  probabilities  P’  and  /]' .  Even  though  Pf 
is  very  small,  the  conditional  probabilities  can  be 
sufficiently  large  to  be  evaluated  accurately  by  MCS 
using  a  much  smaller  number  of  samples  for  the 
intermediate  failure  events  F1  ,  i  =  2,...,m.  For 
example,  if  the  probability  of  failure  for  the  rare  event 
F  is  /]  =1(T6,  the  failure  domains  F'  ,  1  =  2,... ,6  can 
be  defined  so  that  the  corresponding  failure  events  have 
probabilities  /]'  =  0. 1  and  conditional  probabilities 

6 

Pf  =0.1  for  /  =  2,... ,6  such  that  Pf  =  =  1CT6 . 
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The  much  larger  probabilities  P'  can  be  efficiently 
evaluated  using  MCS  for  a  given  accuracy.  The 
probability  calculation  of  a  rare  event  is  therefore, 
replaced  by  a  sequence  of  more  frequent  events  in  the 
conditional  probability  spaces. 

To  evaluate  the  conditional  probabilities  in  Eq. 
(11),  a  sampling  scheme  is  needed  to  produce  sample 
functions,  which  satisfy  the  conditional  events.  Two 
such  sampling  schemes  are  Markov  Chain  Monte  Carlo 
(MCMC)  and  splitting.  Their  difference  is  on  how  they 
generate  offsprings  from  appropriate  parent  samples.  A 
brief  introduction  to  MCMC  and  splitting  is  provided 
next. 

3.1  Sampling  with  MCMC  and  splitting 

Subset  simulation  requires  the  conditional 
probabilities  P'  =  |  for  i  =  1,2,. .  ,,m  to 

calculate  the  overall  small  probability  of  failure  Pt  as 
the  product  of  the  larger  conditional  probabilities 
P'  and  P1 .  Monte  Carlo  simulation  can  estimate  the 
conditional  probabilities  as  long  as  we  can  generate 
conditional  samples.  For  example,  to  calculate 

we  must  generate  sample  functions 
g(x(r),r)  which  satisfy  the  failure  condition 

F'_1  =  <  maxg(X(r),r)>  S','-1 1  .  For  that,  Markov  Chain 

(<e[o,r]  J 

Monte  Carlo  (MCMC)  or  splitting  (S)  sampling 
schemes  can  be  used. 

To  generate  MCMC  sample  functions,  we  perturb  a 
“parent”  sample  function  gp(x(t\t )  which  satisfies  the 

condition  of  the  F1  1  domain  (thick  line  in  Figure  2)  to 
generate  an  offspring  g0(x(r),f)  (thin  lines  in  Figure 
2).  However,  some  of  the  offsprings  may  not  belong  to 
the  F'  1  domain  (thin  solid  line  in  Figure  2)  and  are 
therefore,  rejected.  As  a  result,  the  MCMC  sampling 
process  can  be  inefficient.  More  details  are  provided  in 
[26,30,31], 

Generating  offsprings  using  splitting  does  not  have 
the  drawback  of  MCMC.  Figure  3(a)  highlights  the 
main  idea.  Considering  that  the  “parent”  sample 
function  gp(x(r),r)  belongs  to  the  F'  1  domain,  we 
can  partition  it  into  two  partial  signals 

gp  =  Up (X(t ),?),/<  fpass,  ]  and  g+  =  [gp(x(f),f), 
t>t  1 1,  where  t  ,  is  the  time  when  ep(X(r),n 
crosses  the  threshold  S'  1  and  t  =  t  ,  —  At  is  the 

1  pass  pass 

time  instant  right  before  the  first  passage  time  tpass,  (see 

Figure  3a),  and  record  the  system  states 

S  ,  =  X/  ,  —  At).  We  assume  1)  that  the  sample 

pass  \  pass  /  7  r 


functions  of  the  excitation  random  processes  U(t)  can 
be  partitioned  similarly  into  Up  =  [U(V),  t  <  tp^  ,  ]  and 
Up  =[U(r),r  >fpassi]  and  2)  that  the  new  sample 
functions  Uq  =  [U(r),  t  >  fp|  ^  ]  of  the  excitation  can  be 
generated  after?  ,  .  Then,  using  the  states  S  and 

°  pass  c’  pass 

the  excitation  Uq  ,  we  generate  an  offspring  (gP,go), 
where 

£o  =  U  (x(f )’  Uo  (?),?)  f  >  fpass,  ]  • 

This  process  guarantees  that  the  offspring  (gP,go) 
belongs  to  the  F'  1  domain. 


Figure  2.  Accepted  and  rejected  sample  functions 
generated  by  MCMC 


?(x(dO=s, _ 

So  (Always  accepted) 


S  p 

-► 

fpass  i  t 


g(x(,U)=s;-‘  u>! 


(b) 

Figure  3.  Sample  functions  generated  by  splitting 
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If  the  time  step  At  is  small  enough,  the  difference 
among  pass), ,  pass], ,  and  pass'  is  very  small  and  the 
splitting  point  pass'  (  Figure  3a)  can  be  replaced 

with  pass p  =  pass1  (Figures  3a  and  3b).  Then,  the 
generated  offspring  is  always  accepted  without  a  path 
dependency  issue.  However  if  the  time  step  At  is  not 

small  enough,  pass], ,  passj, ,  and/or  pass 1  may  be  quite 
off  from  the  threshold  ,S't'  and  the  offspring  will  have 

a  path  dependency  issue  if  the  splitting  starts  at  pass],  . 
To  avoid  path  dependency  and  improve  accuracy,  we 
perform  splitting  at  pass1-  similarly  to  [32],  This  may 

compromise  efficiency  because  of  a  repeated  generation 
of  response  realizations  until  the  generated  offspring 
path  has  a  first  passage  point  in  the  time  period  of 
interest. 

4.  PROPOSED  APPROACH 

Subset  simulation  with  splitting  (SS/S)  or  subset 
simulation  with  MCMC  (SS/MCMC),  are  promising 
approaches  in  estimating  the  reliability  of  dynamic 
systems  with  low  probability  of  failure.  Even  though 
SS/S  is  more  efficient  than  SS/MCMC,  its 
computational  efficiency  decreases  rapidly  when 
several  probabilities  of  failure  (e.g.,  Pf  (0,7]) ,  Pf  (0. T, ) 

and  Pf  (0, 7’, )  with  7]  <T 2<  71 )  are  needed  to  build  the 

CDF,  because  each  probability  of  failure  is  estimated 
independently  starting  from  t  =  0  . 

To  address  this  issue  for  non-repairable  systems, 
we  propose  a  subset  simulation  method  with  splitting 
and  partitioned  time  (SS/SPT).  The  method  partitions 
[0,7’]  into  several  time  intervals  and  uses  a  modified 
SS/S  approach  sequentially  using  the  survived  (not 
failed)  sample  trajectories  from  the  previous  time 
interval  to  directly  obtain  the  failure  rate  in  the  current 
time  interval,  and  then  calculate  the  cumulative 
probability  of  failure.  A  nested  sequence  of  m  failure 
regions  Fl  3  F2  zd  ■  ■  ■  zd  F"'  =  F  is  formed  for  different 
threshold  levels  5*  <  S’,2  <  •  •  •  <  S’"  =  St .  The  random 
processes  within  different  time  intervals  are  correlated 
because  if  a  failure  occurs  outside  a  time  interval,  the 
failed  system  is  disposed  at  that  time. 

The  time  of  interest  [0,r]  is  discretized  at 
instances  ]o  71  7)  •••  T  ■■■  T,,  =t\  so  that 

v  1  ^  "  '’interval  7 

N- 

T  -  zk  -r„_1)=  /Vinterva AT  with  r0  =  o,  tNm  =  r 

n= 1 

and  at  =  T/n  ■  A  series  of  7Vintenral 

/  v  interval 

correlated/conditional  failure  regions  are  then  formed 
defined  by  the  following  events 


Fc„  =  |  max  g(x(r),r)>  St  Fn_{  \  , 

(12a) 

where 

Fn-i=\  max  g(x(r),r)<  St  i , 

(12b) 

and 

Fn=\  max g (x(f ), r) >  A, } 

(12c) 

for  «  =l,2,...,A/interval. 

The  dependencies  in  the  conditional  failure  events 
Fcn  in  the  time  direction  are  simply  enforced  during 
simulation  by  using  only  the  surviving  (not  failed) 
sample  functions  among  all  sample  functions  in  [0,7’n_1] 
in  the  subsequent  time  period  \Tn  , ,  7]  ] .  This  allows  us 
to  partition  time  T  in  shorter  time  intervals 
[T„-vT„  ]•  n  =  1,2, . . .,  A/interval  and  use  the  subset 
simulation  with  splitting  in  each  of  the  time  intervals  to 
calculate  the  conditional  probabilities 

Pcfn=P(Fcn)-  (13) 

For  the  time  interval  [?j  ,,/j],  we  use  F‘n  instead 

of  F1  to  represent  the  /th  subset  event  in  the  threshold 
direction 

Fl=\  max  g(x(r),r)>S;|/U,  (14) 

where  i  =  1,2,. .  ,,m„ ,  and  mn  is  the  number  of 
intermediate  thresholds  for  the  nrh  time  interval 
[/)  ,,/j].  The  corresponding  probability  is  denoted 

The  probability  Pcfn  =  f(fc„)  is 

mn 

calculated  using  Eq.  (11),  where  Fcn  =  P|Fc'n  . 

1=1 

The  probability  Pcfn  can  be  also  used  to  calculate 
the  failure  rate 


in  the  time  period  [?j  ,,7j]  where  ATn=Tn—Tn_l, 
because  it  represents  the  probability  of  failure  in  the 
time  period  [?)  ,,7)]  under  the  condition  that  the 
system  is  safe  in  [0,  Tn_{  ] . 

Because  Fn  =  Fn  l  U we  can  use  the 
modified  SS/S  approach  to  calculate  Pctn  in 
[7j  i,7j]  and  the  probabilities  of  failure 

Ptn=P(Fn)=Ptn-1+(l-Pfnhn  d6) 
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at  all  instances  Tn  ,  n  =  l,2,...,7jnterval  with  almost  the 
same  computational  effort  the  SS/S  approach  requires 
to  calculate  the  probability  of  failure  at  the  final  time  T. 
Details  on  the  computational  effort  are  provided  in 
Section  4. 1 . 

Figure  4  highlights  the  basic  idea  of  the  SS/SPT 
method.  At  most  N  trajectories  are  simulated  in  [0,r]. 
The  simulation  is  first  carried  out  in  [0,7j].  For 
illustration  purposes.  Figure  4  shows  only  two 
trajectories.  Because  both  trajectories  cross  the 
intermediate  threshold  Sj1,  they  are  used  to  generate 
offsprings  (dotted  lines)  in  [0,7j].  The  offsprings  along 
with  the  parent  trajectories,  are  used  to  calculate  the 
probability  of  crossing  the  next  intermediate  threshold 

Sj2 .  The  system  states  X(r)  at  points  A  and  B  are 
recorded  in  order  to  continue  if  necessary,  the 
simulation  of  the  particular  trajectory  for  the  next  time 
interval  [7j,r2],  The  process  is  repeated  for  all 
intermediate  thresholds  in  [0,7j].  All  offsprings  are 
generated  only  within  the  current  time  interval. 

Because  trajectory  1  crosses  the  overall  threshold 
Sj  in  [0,7j],  it  is  considered  failed  and  removed  from 
the  population.  Thus,  we  do  not  continue  its  simulation 
in  [7j,T2] .  In  contrast,  because  trajectory  2  did  not  cross 

St  in  [0,7j],  we  continue  its  simulation  in  [7j ,  7'2  ]  using 
the  states  X(f)  at  point  B.  In  [7j,7j],  the  second 
trajectory  crosses  both  S,1  and  S2,  and  can  produce 
offsprings  as  shown  in  Figure  4.  The  process  of 
continuing  the  simulation  of  only  the  survived  parent 
trajectories,  which  can  subsequently  create  offsprings, 
is  carried  out  for  all  time  intervals. 


Figure  4.  Schematic  of  basic  features  of  the  SS/SPT 
method 

A  step-by-step  description  of  the  proposed  SS/SPT 
algorithm  is  given  below  followed  by  the  developed 
SS/SPT  algorithm. 

Step  1.  Initialize  algorithmic  parameters. 

Step  2.  Generate  N  =  Nsn  replications  of  the 
system  random  parameters. 

For  each  time  interval.  Steps  3  through  8  are 
performed: 


Step  3.  Generate  a  state  trajectory  for  each  system 
parameter  replication  for  a  random  excitation 
trajectory. 

Step  4.  Determine  the  first  threshold  level  for  failure 
event  Fl  (see  Equation  12c),  and  calculate  the 
corresponding  probability  of  failure  Pf  =  F^F1). 

Step  5.  Use  splitting  to  generate  offspring  trajectories 
from  randomly  selected  “parent"  trajectories  for  failure 
event  F' . 

Step  6.  Determine  the  first  threshold  level  of  failure 
event  F1+‘  and  calculate  the  corresponding  probability 
of  failure  Pff‘  =  f(f1+1). 

Step  7.  Calculate  the  failure  rate  and  cumulative 
probability  of  failure. 

Step  8.  Reset  for  next  time  interval. 


A  pseudo  code  description  of  the  proposed  SS/SPT 
algorithm  is  provided  below  with  details  for  all  the 
above  steps. 

Step  1  :  Initialize  algorithmic  parameters. 

Specify  the  guessed  probability  of  failure  /jj’""' 

( If,  in  Equation  16). 

Specify  the  guessed  expected  conditional 
probability  pf“  "  for  each  subset  simulation  level. 

Determine  the  expected  conditional  probability  p 
for  each  subset  simulation  level 

p  =  max{prss,Pfn“ess}. 

Specify  the  number  of  the  initial  independent 
parent  response  trajectories  (sample  functions)  Nsn  . 

N  =  Nm; 

Specify  the  life  span  of  interest  T  . 

Set  T0  =  0  . 

Specify  a  uniform  time  interval  7jntervalto  partition 

T . 


Calculate 


Nt„ 


T  =  N: 


intervar  interval* 


and  reset 


Specify  the  initial  values  of  excitation  for 

k=l2,...JV^. 


Specify  the  initial  values  of  the  states  X^  for 

*=1,2,...^. 


Set  Pf0  =  0  by  assuming  that  the  system  has  a  very 
small  failure  probability  at  the 
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initial  time  T0  =  0  . 

Step  2:  Generate  N  =  Nm  replications  of  the  system 
random  parameters. 

For  k  =  1 :  IV 

Generate  independent  random  system 
parameters  . 

End 

For  n  =  1 :  /V,nterval 

Step  3  :  Generate  one  state  trajectory  for  each 
system  parameter  replication  for  a  random 
excitation  trajectory. 

For  k  =  1 :  N 

For  time  period  [7^,7],] generate 
independent  random  excitation  parameters 
Ye<A)  and  the  corresponding  excitation 

trajectory  U(A)  (f ) ;  set  Ylk>  =  [y^  ;  Yew  ] ; 
solve  for  state  trajectory  (sample  function) 
X(k)(t)  using  Equation  (4)  with  initial 
conditions  and  x[A| .  Calculate  the 

corresponding  response  g(xa’(f),r) . 

End 

Let  Ns  be  the  number  of  trajectories  that 

satisfy  j  max  ga>(x(r),r)<  St,  k  =  1,2,. ..,7^1  . 

Denote  and  record  the  corresponding 
excitation,  state,  and  system  parameter  vectors 

by U^(?)  ,  X^‘*(f)  and  respectively, 

for  ks  =1,2,...,NS . 

Step  4  :  Determine  the  first  threshold  level  for  failure 
event  F 1  (see  Eq.  12c),  and  calculate  the  corresponding 
probability  of  failure  Pf'n  = 

i  =  0; 

Let  S,1^  be  the  (l-p)xlOO  percentile  of 

max  g<A'(x(r),r),  k  =  1,2,..., iV 

tfPn-1  ,r„] 

Record  the  time  coordinates  and  state  values 
of  the  trajectories  (sample  functions)  at  the 
first-passage  points  over  the  intermediate 
threshold  g(x<k>(t),t)-  S'f  and  record  the 
corresponding  system  parameters.  Instead  of 
using  the  approximation  Nf  ~  pN  , 
calculate  N{  by  counting  the  exact  number  of 

all  trajectories  satisfying  g(fCk>(t),t)>  1S’t1+' . 
Denote  the  time  coordinates,  state  values,  and 


the  corresponding  system  parameters  right 
before  the  first  passage  times  by 

1ffcl,Sfe)1.,Yft)1  \k{  =U,..JVf|  and 

ll  pass  pass  spass  J  1  1  ) 

record  the  trajectories  jx^,  kt  =  l,2,...,iVf  ] 
reaching  the  threshold. 


Calculate  Pt'jj  =  ^ < 


N ’ 


While  .S’,1'  <  .S’, 
i=i+1; 

Step  5  :  Use  splitting  to  generate  offspring 
trajectories  from  randomly  selected  “parent” 
trajectories  for  failure  event  F‘ . 

For  k  =  (Nf  +l):N 

Randomly  select  values  for  the  tuple 

[cS.yJ 

from 

,  S(*f), ,Y<*f) ,  J  kf  =  1,2,... Wf) 

<L  pass  pass  spass  J  1  1  J 

with  a  uniform  probability  of  l/Nf  . 


X<k>(t),  t  =  [t,t  +Ar]} 


Generate 

starting  from  S  using  Ys  and  the 
newly  generated  excitation  U(t) 
according  to  Equation  (4). 
If  g(x(f  +  At), t  +  At) <  S'  ‘ ,  redo  this 
step  for  the  same  k 
until  ^(x(f  +  Af),f+A?)>  S, 1+1 .  For 
the  first  passage,  denote  the  time 
coordinate,  state  value,  and  the 
corresponding  system  parameters  by 

|fW,,SW„Yw  ,  I  (noting  that 

L  pass  pass  spass  J  ° 

Ak) 


r  \  =  t  +  At  and  Y(  } ,  =  Y„ ).  Record 

pass  pass  s  ' 

this  one  time  step  trajectory  as  X^A  * . 
Continue  the  trajectory  X*A>  by 
generating 


pass 
( k ) 


starting  from  S(  ,  using  Yu,,and 

°  pass  °  pass 

the  newly  generated  excitation  U!<  l(r) 
according  to  Equation  (4). 

End 

Step  6:  Determine  the  first  threshold  level  of 
failure  event  F1+l  and  calculate  the  corresponding 
probability  of  failure  P)f  =  f(f1+'). 

Let  .S’,1  '  he  the  (l-p)xlOO  percentile  of 


max  ga,('X.(t\t\  k  =  1,2,...,N 

Record  the  system  parameters,  time 
coordinates  and  state  values  of  the 
trajectories  exactly  at  the  first-passage 
over  the  intermediate 

threshold  g(x(A,(t),f)=  St1+'  .  Instead  of 

using  the  approximation  Ns  ~  pN  , 
calculate  Nf  by  counting  the  exact  number 
of  all  trajectories  satisfying 
g(xa>(t),tj  >  S'*‘ .  Denote  the  first 
passage  time  coordinates,  state  values  and 
the  corresponding  system  parameters  by 
L(k, )  s(k, )  y(Af \\kf  =  1,2, . . . and 

record  the  trajectories 

|x^f  \kf  =  1,2, . . .  JV{ }  reaching  the 

threshold. 

Calculate  /£'=% 

End 
i  =  i+ 1 

Step  7_:  Calculate  the  failure  rate  and  cumulative 
probability  of  failure. 

Let  Nf  be  the  number  of  trajectories  that 
satisfy 

max  gm(x(t),t)>S„k=l,2,...JV 

>4Tn_,,T„] 


m 

i= 1 


Step  8:  Reset  for  next  time  inten’al. 

N  =  NS ; 

U?)  =  Uj(rJ,xW=X«(rB)  and  YSW  =  YSW 
for  A;  =  1,2,..., ./V; 

Revise  the  expected  conditional  probability  p 
for  each  subset  simulation  level 

p  =  max{pr\Pr} 

End 


4.1  Computational  Effort 

At  the  n th  time  interval,  we  simulate  Nsn  parent 
trajectories  within  [jj  , , 7j; ]  by  solving  the  dynamic 
Equations  (4).  Using  a  constant  time  integration 
step  At ,  the  number  of  function  evaluations  (see 
paragraph  above  Equation  5  for  definition)  per 

~T„ t  T 

&  ^,„tervaAf 

where  YniIcrv.i]  is  the  number  of  time  intervals  we 
partition  the  total  time  [0,r]  into.  This  results  in  an 
T 

upper  limit  Ns  -  of  function  evaluations  for 

^intervaA  # 

the  parent  trajectories. 

At  the  i'h  threshold  of  the  nth  time  interval,  we 
simulate 

-NL)=(Nsn-pnNsn)={l-pn)Nsn  (17) 
offspring  trajectories  in  \l'n  ,,7’J  from  t  to  Tn  by 
solving  the  dynamic  Equations  (4)  where  t  e  \_I'H  ,,7j]  is 
the  time  a  trajectory  crosses  the  ;th  threshold.  In 
Equation  (17),  N‘tn  is  the  number  of  the  sample 


trajectory  within  [c  , , Tn  is 


functions  among  Nsn  which  cross  the  ;th  threshold  and 
pn  is  the  desired  conditional  probability  for  each 
subset  level.  Based  on  the  definition  of  t  ,  it  is 
T  +T 

reasonable  to  assume  t  ~  — — - —  .  The  upper  limit  of 


function  evaluations  per  offspring  trajectory  at  the  ; 
threshold  of  the  nth  time  interval  is  therefore,  equal  to 
T-t  T-T  ,  T 


At  2At  2AintervaA/ 


- .  If  m.  is  the  number  of 


thresholds  where  offsprings  are  generated,  we  need  up 


to  —  —  (mn—l)  function  evaluations 

2VntervaAf 

for  all  (m„-l)  thresholds. 

The  upper  limit  of  the  total  number  of  function 
evaluations  for  the  nth  time  interval  [yj  ,,7’J  is  thus 


equal  to  Nsn 


WintervaA / 


-(i-zaK,, 


2Vn,ervaA  t 


,  ,  (1-P„) K-l) 
2 


N 


-At... 


(18) 


1  intervaA^ 

The  number  of  function  evaluations  in  relation  (18) 
is  approximately  the  same  for  all  A(lllcivaltime  intervals. 


yielding  the  following  upper  bound  for  the  overall 
number  of  function  evaluations 
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N.  + 

1  v  intervall  A  '  2 


^n.ervaA' 


-N" 


x  ,  ^-Pn\mn  -l) )  T 


(19) 


—  A? 
Ar 


Considering  that  represents  the  number  of 


function  evaluations  to  generate  a  single  trajectory 
(sample  function)  of  the  output  random  process  from  0 
to  T,  the  upper  limit  of  an  equivalent  number  of 
simulated  sample  functions  is 


(20) 


Relation  (18)  provides  the  approximate  number  of 
function  evaluations  for  the  subset  simulation  with 
splitting  (SS/S)  method  if  ATinterval  =  1 .  This  results  in 
the  same  computational  effort  with  the  proposed 
SS/SPT  approach  (Equation  20)  in  terms  of  equivalent 
number  of  simulated  sample  functions.  However,  the 
SS/SPT  method  calculates  the  entire  CDF  and  the 
failure  rate  as  a  function  of  time  while  the  SS/S  method 
simply  provides  the  probability  of  failure  at  the  final 
time  T.  This  can  be  a  significant  advantage  of  the 
proposed  SS/SPT  over  the  direct  SS/S. 


4.2  Estimation  of  Number  of  Parent  Trajectories 
and  Intermediate  Thresholds 

This  section  estimates  the  bounds  of  the  coefficient 
of  variation  (C.O.V.)  8P  of  the  probability  estimate 

P.j  ri  of  Pcin  (see  Equations  13  and  12)  for  the  time 
interval  [Tn_1,Tn],n  =  l,2,...,Nints;r,al  where T0  =  0  .  The 
bounds  are  used  to  estimate  the  appropriate  number  of 
parent  trajectories  Nsn ,  the  expected  conditional 
probability  of  failure  pn  for  each  subset  level,  and  the 
number  of  intermediate  thresholds  mn  for  each  time 
interval  for  a  given  C.O.V. 

In  subset  simulation,  the  C.O.V.  8P  of  P ,  is 

rcfn  CIn 

approximated  by 

mn 

(^J2*Z(<V)2’  (21) 

i= 1  Cf” 

where  Spi  is  the  C.O.V.  of  the  probability  estimate 

°cf  n 


P^n  at  the  threshold.  According  to  [31],  8  ,  is 


given  as 


S*.  = 


1  />: 


(22) 


In  our  proposed  SS/SPT  approach,  yln  =  0  and  y‘n 
is  a  constant  provided  by 

r 

N' 


o  <y'n<E\ 


N‘  -  N’r 


(23) 


for  i  =  2,...,mn,  where  N'„  is  the  total  number  of 

sample  functions  for  failure  domain  F'  1  and  N'tn  is 

the  number  of  sample  functions  among  N'n  which  are 

also  used  in  the  failure  domain  F'  .  For  simplicity,  we 
assume  that  each  subset  level  has  the  same  expected 

conditional  probability  Pcjn  =  pn  .  We  also  assume  that 


the  total  number  of  sample  functions  N'n  is  the  same 


for  each  subset  level  i  and  equal  to  Nsn  (i.e.,  N'n  =  Nsn ) 

which  must  be  estimated.  This  assumption  results  in  the 
same  number  of  conditional  failure  samples 

Nf„  =  P„  A)  „  for  each  subset  level.  Therefore,  the 
expected  number  of  offspring  trajectories,  generated 
from  the  N'fn  sample  functions,  and  reaching 


the  i  threshold  is 

NL=Nsn-Nltn={l-p„)Nsn, 
and  Equation  (23)  becomes 

N 

1 


0  <y'„<E 


'  l  ' 

1  -  Pn 


Pn 


(24) 

(25) 


for  i  =2 ,...jnn 

Combining  Equations  (21),  (22)  and  (25)  and  using 
the  assumption  P’fn  =  p„ ,  the  C.O.V.  Sp  of  the 


probability  estimate  Pcfn  for  the  n"  time  interval,  has 
the  following  bounds 


X(^)2<(^,„)2<(^)2+Z(^): 

i= 1  i=2 

where 


1  +  - 


1 


l~Pn 


3.  = 


1 ZPn 

PnNSn 


,  (26) 


(27) 


Substitution  of  Equation  (27)  in  Equation  (26) 
yields 

r  “i 

1  ~Pn 


P„NS„ 


fej2 


1  ~Pn 


P„N,n 


(28) 


Based  on  Equation  (11)  and  the  assumption 
P‘fn  =  pn ,  the  number  of  threshold  levels  mn  for  the  n,h 
time  interval  in  Eq.  (28)  can  be  calculated  from 

mn  mn 

Potn=I\PL=Y\Pn={PnY 

i= 1  i= 1 


or 


ln(/c,  J  =  mn  lr|(/V)  =>  in,,  =  ceil 


(29) 


V  lnPn  J 

Equations  (28)  and  (29)  are  used  to  estimate  the 
bounds  of  8V  based  on  estimates  of  Ns„  and  Pcfn ,  and 
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a  predetermined  value  of  pn .  Alternatively  for  a  given 
8Pi  ,  we  can  estimate  the  bounds  of  Nsn . 

As  an  example,  we  assume  that  the  time  internal  of 
interest  is  [o,r]  =  [o,  352  sec]  and  the  estimated 

probability  of  failure  Pfn  =  P(Fn)  (see  Eq.  12c)  at  T  = 


352  sec  is  Pfn  (0,352)  =  0.2  .  If  we  partition  the  interval 
[o,r]=  [o,352sec]  at  every  8  seconds  (i.e., 
352 

-^interval = - =  44  time  intervals),  Pcfn  is  estimated  by 


P,  =  - 


1  fn 


0.2 


N-  ,  44 

v interval 


=  0.004545  .  If  we  also  assume 


pn  =  0.2 ,  Equation  (29)  estimates  the  number  of  subset 


levels  m„  =ceil 


"'H  V 

lnP„ 


=  ceil 


ln(0.004545) 

ln(0.2) 


=  4  for 


each  time  interval.  Then  Equation  (28)  can  provide  the 
bounds  of  /V,  „  for  an  assumed  C.O.V.  of 


5P  =0.07  as 


[~P„ 

pJSpJ2 


<  N„,  < 


4  ^  <  N  < 

0.2  -0.07 2 


l  +  (m„-l)2  P 
l  +  (4-l) 


1  -P„ 
2-0.2 


1  ~Pn 


P-VpJ 

1-0.2  _ 


0.2 -0.07 2 


1-0.2 

3265  <  /VS7!  <  6326 .  Choosing  Nm  =  5000 ,  the  relation 
(20)  gives  the  approximate  number  of 


AL  =  1  + 


i !  Q-p„Xw„  -Q 

2 

=  1 1,000  simulated  trajectories. 


(1-0.2)  -(4-1) 


5000 


5.  VEHICLE  VIBRATION  EXAMPLE 

The  linear  representation  of  a  quarter  vehicle  in 
Fig.  5  is  employed  to  illustrate  the  proposed  SS/SPT 
method.  The  vehicle  travels  over  a  stochastic  terrain  at 
20  miles  per  hour  (mph).  The  system  has  two  random 
parameters;  the  damping  coefficient  bs  and  the 

stiffness  coefficient  ks .  Both  are  normally  distributed 
with  bs  ~  A(7000,14002)N/m/s  and  ks  ~ 
A''(40000.40002)  N/m.  The  sprung  and  unsprung 
masses  ms  and  mu  of  the  suspension,  the  tire  stiffness 
coefficient  kt  and  the  tire  damping  coefficient  b]  are 
deterministic  with  ms  =  100kg  ,  ma  =  100kg  ,  kt  =  40 
x  104N/m  and  kt  =  40xl04N/m,  bt  =40xl03 
N/m/s .  The  vehicle  is  travelling  over  a  stochastic 
terrain  which  provides  excitation.  The  elevation  n(r)of 
the  terrain  is  a  random  process. 

The  state-space  approach  is  used  to  determine  the 
vertical  acceleration  g(x(r), t)  =  xs(t)  of  the  sprung 
mass,  in  g’s.  Failure  occurs  if  the  magnitude  of  the 


vertical  acceleration  exceeds  a  threshold  St  (i.e., 

\g{x(t),t)>St). 


Figure  5.  Quarter  vehicle  model 

A  third  order,  autoregressive  time  series  model 
AR(3)  is  used  to  characterize  the  road  excitation 
process  [25].  Time  series  modeling  can  characterize 
both  stationary  and  non-stationary  processes  [34,  35], 
The  road  model  for  this  example  is  expressed  as 
U  -  =  1  ,2456mh  - 0.2976m ;_2  -0.1954m ._3  +  s.(o,0.51322) 
where  the  subscript  j  indicates  the  time  step  of  the 
discretized  time,  f/(o,051322)  is  Gaussian  white  noise 

with  a  standard  deviation  of  0.5132.  The  coefficients 
1.2456,  -0.2976,  and  -0.1954  are  the  three  estimated 
feedback  parameters,  of  the  AR(3)  model. 

The  following  two  differential  equations  comprise 
the  equations  of  motion 

mX  +  ibt  +  h-.  k  -  b,K  +  k  +  K  k  -  Kxs  =  Ku + hP 
+  b  ks  -  k )+ K  ks  -  ■ ) = 0 

which  are  transformed  to  state-space  form  and 
integrated  in  time  to  obtain  the  response  ( x.  (t)  and 
x[A  (t))  which  is  then  used  to  determine  the  vertical 
acceleration  xs(f )  of  the  sprung  mass. 

The  time  step  for  the  numerical  integration  of  the 
above  equations  is  equal  to  A?  =0.01  seconds.  The 
threshold  for  failure  is  St  =3.5  g’s.  We  calculated  the 
failure  rate  and  the  CDF  in  the  time  period  [0, 352] 
seconds.  The  latter  was  partitioned  every  2  seconds 
( interval  = 176  X  4  seconds  ( TV interval  =  88 )  and  8  seconds 
( /Vlntcrvai  =  44 ).  For  the  2  second  partitioning  for 
example,  the  probability  of  failure  is  calculated  at  2,  4, 
8,  ...  ,  352  seconds. 

Following  the  process  in  the  last  paragraph  of 
Section  4.2,  we  can  estimate  the  number  of  parent 


11 


trajectories  Nsn  ,  the  number  of  intermediate  thresholds 
(or  subset  levels)  mn  and  an  approximate  number  of 
simulated  trajectories  for  a  desired  C.O.V.  of 
SP  =0.07  and  expected  conditional  probability 

pn  =  0.2 .  If  the  time  [0,  352]  seconds  is  partitioned 
with  an  8  sec  interval,  we  obtain  408 l  <Nsn  <8163, 
mn  =  5  and  an  approximate  number  of  15,000  simulated 
trajectories.  For  partitioning  with  4  or  2  seconds,  we  get 
3265  <  Ns  n  <  6326  ,  mn  =  4  and  an  approximate 
number  of  1 1,000  simulated  trajectories.  Based  on  these 
estimates,  we  selected  Ns  n  =  5,000 . 

Figure  6  and  Table  1  show  the  calculated  failure 
rate  and  CDF  with  a  step  of  2  seconds.  Four  different 
runs  were  performed  (red  dotted  lines)  in  order  to 
demonstrate  the  expected  variability.  Table  1  provides 
detailed  results  for  the  four  independent  runs,  their 
mean  and  the  C.O.V.  It  also  shows  the  MCS  estimates. 
The  black  dotted  line  represents  the  average  of  the  four 
runs.  The  estimated  probabilities  from  MCS  with  106 
trajectories  (blue  solid  line)  are  used  as  a  baseline. 
Figure  6  shows  that  the  average  failure  rate  and  CDF  of 
the  four  SS/SPT  runs  is  close  to  the  MCS  estimates. 
Figures  7  and  8,  and  the  corresponding  Tables  2  and  3, 
show  the  calculated  failure  rate  and  CDF  with  a  step  of 
4  and  8  seconds,  respectively. 

Figures  6a,  7a  and  8a  show  that  the  estimated 
failure  rate  l(t)  oscillates  around  its  mean  throughout 
time.  The  variability  is  due  to  the  low  value  of  A(t)  and 
the  relatively  low  number  of  simulated  trajectories. 
However,  the  estimation  is  unbiased  since  the 
oscillations  are  about  the  mean  value.  The  unbiased 
estimation  of  the  SS/SPT  algorithm  is  also  supported  by 
the  estimated  CDF  from  the  four  runs  in  Figures  6b,  7b 
and  8b  which  show  variation  around  their  mean.  The 
CDFs  do  not  exhibit  an  oscillation  because  the  number 
of  parent  sample  functions  Nsn  was  estimated  based  on 

the  low  C.O.V.  of  0.07  for  each  time  interval.  The 
probability  variation  is  expected  to  be  small  because  the 

increase  in  the  probability  of  failure  Pcfn  at  each  time 
interval  is  small  (approximately  equal  to 
0.2 

— ^  =  0.00 1 1 36  for  the  2  second  interval  case)  and  its 
176 

COV  §P  is  chosen  small.  This  is  an  advantage  of  the 

proposed  SS/SPT  approach  over  the  direct  SS/S 
approach  where  the  estimated  CDF  can  exhibit  large 
oscillations. 

Figure  9  shows  the  COV  of  the  estimated  failure 
rate  A(t)  and  probability  of  failure  Pf .  We  observe 

that  the  COV  of  A (f)  decreases  on  average,  for 
increasing  time  interval  (0.3364  for  2  sec,  to  0.2459  for 
4  sec,  to  0.1876  for  8  sec).  This  is  because  the  error  of 
probability  of  failure  estimation  for  any  MCS -based 


method  reduces  as  the  probability  of  failure  increases. 
As  we  increase  the  time  interval  in  the  proposed 
SS/SPT  method,  the  probability  of  failure  contribution 
from  the  larger  time  interval  to  the  overall  P, 

increases,  reducing  therefore  the  estimation  error. 

The  same  observation  holds  for  the  C.O.V.  of 
CDF.  Figure  9b  shows  that  the  C.O.V.  decreases  on 
average,  for  increasing  time  interval  (0.048  for  2  sec,  to 
0.0314  for  4  sec,  to  0.0286  for  8  sec).  The  average 
C.O.V.  of  the  CDF  is  much  smaller  than  the  C.O.V.  of 
A  (?)  due  to  the  higher  value  of  Pf  compared  to  the 

value  of  A(t) .  Furthermore,  the  C.O.V.  of  the  CDF 
decreases  rapidly  in  the  early  period  and  exhibits  less 
variation  at  later  times  because  the  CDF  value  increases 
with  increasing  time. 

Based  on  the  above  discussion,  a  relatively  large 
time  interval  (e.g.  8  versus  2  seconds  for  this  example) 
can  be  used  in  order  to  reduce  the  failure  rate  and  CDF 
variation. 


t  (secs) 
(a) 


Figure  6.  Failure  rate  (a)  and  probability  of  failure 
(b)  from  SS/SPT  method  calculated  at  2  sec  intervals 
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Failure  rate:  h  Failure  rate: 


(a) 


(b) 

gure  7.  Failure  rate  (a)  and  probability  of  failure  (b) 
from  SS/SPT  method  calculated  at  4  sec  intervals 


4.5 
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Figure  8.  Failure  rate  (a)  and  probability  of  failure  (b) 
from  SS/SPT  method  calculated  at  8  sec  intervals 


(a) 


(b) 

Figure  9.  C.O.  V.  of  failure  rate  (a)  and  probability  of 
failure  (b)  from  SS/SPT  method  for  2,  4,  and  8  sec 
intervals 
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Table  la.  Details  of  from  SS/SPT  method  calculated  at  2  sec  intervals 


Time 

(sec) 

SS/SPT 

MCS 

Run  1 

Run  2 

Run  3 

Run  4 

Mean 

c.o.v. 

2 

6.93E-04 

7.04E-04 

8.47E-04 

8.84E-04 

7.82E-04 

0.124953 

8.38E-04 

4 

8.24E-04 

1 .20E-03 

1.02E-03 

9.51E-04 

1.00E-03 

0.158611 

9.12E-04 

6 

7.83E-04 

4.62E-04 

9.77E-04 

7.33E-04 

7.39E-04 

0.287182 

8.93E-04 

8 

9.38E-04 

6.69E-04 

7.49E-04 

7.56E-04 

7.78E-04 

0.146443 

8.79E-04 

10 

1.06E-03 

9.83E-04 

1.15E-03 

5.93E-04 

9.48E-04 

0.260010 

9.24E-04 

12 

9.23E-04 

1.25E-03 

7.96E-04 

7.56E-04 

9.31E-04 

0.240193 

8.53E-04 

14 

8.47E-04 

8.98E-04 

1.50E-03 

8.39E-04 

1.02E-03 

0.315475 

8.77E-04 

338 

6.30E-04 

5.12E-04 

4.45E-04 

9.05E-04 

6.23E-04 

0.325363 

5.10E-04 

340 

5.99E-04 

2.88E-04 

5.33E-04 

4.83E-04 

4.76E-04 

0.282164 

4.73E-04 

342 

4.37E-04 

4.68E-04 

8.69E-04 

4.35E-04 

5.52E-04 

0.383182 

4.87E-04 

344 

2.78E-04 

1.14E-04 

3.31E-04 

2.62E-04 

2.46E-04 

0.378173 

4.94E-04 

346 

4.30E-04 

3.76E-04 

4.81E-04 

9.16E-04 

5.51E-04 

0.449211 

4.69E-04 

348 

7.50E-04 

2.81E-04 

4.18E-04 

4.08E-04 

4.64E-04 

0.431984 

4.81E-04 

350 

5.53E-04 

5.66E-04 

8.22E-04 

7.32E-04 

6.68E-04 

0.196244 

5.05E-04 

Table  lb.  Details  of  PF  from  SS/SPT  method  calculated  at  2  sec  intervals 


Time 

(sec) 

SS/SPT 

MCS 

Run  1 

Run  2 

Run  3 

Run  4 

Mean 

C.O.V. 

2 

0.001385 

0.001408 

0.001694 

0.001767 

0.001564 

0.124953 

0.001676 

4 

0.003031 

0.003813 

0.00374 

0.003666 

0.003563 

0.100867 

0.003497 

6 

0.004592 

0.004735 

0.005686 

0.005126 

0.005035 

0.097234 

0.005276 

8 

0.00646 

0.006066 

0.007176 

0.00663 

0.006583 

0.069924 

0.007024 

10 

0.008575 

0.008021 

0.009461 

0.007807 

0.008466 

0.087164 

0.008859 

12 

0.010405 

0.010498 

0.011038 

0.009307 

0.010312 

0.070385 

0.010549 

14 

0.012081 

0.012275 

0.014012 

0.010969 

0.012334 

0.101987 

0.012285 

338 

0.191062 

0.204623 

0.188228 

0.205901 

0.197454 

0.046115 

0.197184 

340 

0.192032 

0.205081 

0.189094 

0.206668 

0.198219 

0.045125 

0.197944 

342 

0.192738 

0.205825 

0.190503 

0.207358 

0.199106 

0.043766 

0.198725 

344 

0.193187 

0.206006 

0.191038 

0.207773 

0.199501 

0.043141 

0.199516 

346 

0.19388 

0.206603 

0.191816 

0.209224 

0.200381 

0.043937 

0.200267 

348 

0.195089 

0.207048 

0.192491 

0.209869 

0.201124 

0.042822 

0.201037 

350 

0.195979 

0.207946 

0.193819 

0.211026 

0.202193 

0.042339 

0.201844 

Table  2a.  Details  of  A (f )  from  SS/SPT  method  calculated  at  4  sec  intervals 


Time 

(sec) 

SS/SPT 

MCS 

Run  1 

Run  2 

Run  3 

Run  4 

Mean 

C.O.V. 

4 

9.88E-04 

7.89E-04 

8.75E-04 

1.03E-03 

9.22E-04 

0.119925 

8.74E-04 

8 

6.98E-04 

1 .04E-03 

9.21E-04 

7.50E-04 

8.52E-04 

0.185275 

8.85E-04 

12 

6.85E-04 

1.05E-03 

5.91E-04 

7.01E-04 

7.56E-04 

0.263440 

8.87E-04 

16 

9.68E-04 

8.36E-04 

8.21E-04 

1.02E-03 

9.10E-04 

0.106454 

8.51E-04 

20 

6.34E-04 

8.40E-04 

5.01E-04 

6.70E-04 

6.61E-04 

0.210439 

8.92E-04 

24 

1.22E-03 

6.63E-04 

9.89E-04 

9.58E-04 

9.58E-04 

0.240326 

8.52E-04 

28 

9.44E-04 

6.44E-04 

7.14E-04 

8.43E-04 

7.86E-04 

0.169840 

8.39E-04 

328 

4.48E-04 

5.90E-04 

7.23E-04 

5.65E-04 

5.82E-04 

0.194335 

4.83E-04 

332 

5.62E-04 

5.57E-04 

3.18E-04 

4.26E-04 

4.66E-04 

0.250755 

5.04E-04 

336 

4.58E-04 

4.88E-04 

8.60E-04 

4.30E-04 

5.59E-04 

0.361246 

5.13E-04 

340 

5.66E-04 

4.66E-04 

4.73E-04 

4.44E-04 

4.87E-04 

0.110792 

4.91E-04 

344 

3.39E-04 

7.74E-04 

3.30E-04 

3.79E-04 

4.55E-04 

0.468231 

4.90E-04 

348 

6.32E-04 

6.91E-04 

3.00E-04 

4.54E-04 

5.19E-04 

0.341906 

4.75E-04 

352 

2.79E-04 

3.80E-04 

4.44E-04 

5.41E-04 

4.11E-04 

0.267505 

5.41E-04 
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Table  2b.  Details  of  PF  from  SS/SPT  method  calculated  at  4  sec  intervals 


Time 

(sec) 

SS/SPT 

MCS 

Run  1 

Run  2 

Run  3 

Run  4 

Mean 

c.o.v. 

4 

0.003953 

0.003157 

0.003502 

0.004134 

7.82E-04 

0.119925 

0.003497 

8 

0.006734 

0.007309 

0.007173 

0.00712 

1.00E-03 

0.034798 

0.007024 

12 

0.009457 

0.011459 

0.00952 

0.009904 

7.39E-04 

0.092894 

0.010549 

16 

0.013292 

0.014765 

0.012772 

0.013931 

7.78E-04 

0.062751 

0.013917 

20 

0.015796 

0.018073 

0.014752 

0.016575 

9.48E-04 

0.085815 

0.017437 

24 

0.020617 

0.020676 

0.018651 

0.020342 

9.31E-04 

0.047741 

0.020785 

28 

0.024316 

0.023199 

0.021455 

0.023644 

1.02E-03 

0.052766 

0.024072 

328 

0.199072 

0.200064 

0.198854 

0.191648 

6.23E-04 

0.019640 

0.193087 

332 

0.200872 

0.201847 

0.199874 

0.193026 

4.76E-04 

0.020115 

0.194714 

336 

0.202335 

0.203406 

0.202626 

0.194415 

5.52E-04 

0.020983 

0.196365 

340 

0.20414 

0.20489 

0.204133 

0.195845 

2.46E-04 

0.021192 

0.197944 

344 

0.205219 

0.207351 

0.205185 

0.197063 

5.51E-04 

0.022297 

0.199516 

348 

0.20723 

0.209541 

0.206139 

0.19852 

4.64E-04 

0.023247 

0.201037 

352 

0.208115 

0.210742 

0.20755 

0.200253 

6.68E-04 

0.021751 

0.202766 

Table  3a.  Details  of  A (f )  from  SS/SPT  method  calculated  at  8  sec  intervals 


Time 

(sec) 

SS/SPT 

MCS 

Run  1 

Run  2 

Run  3 

Run  4 

Mean 

C.O.V. 

8 

8.14E-04 

9.55E-04 

7.70E-04 

6.62E-04 

8.00E-04 

0.151619 

8.78E-04 

16 

9.61E-04 

8.30E-04 

8.92E-04 

8.35E-04 

8.79E-04 

0.069234 

8.68E-04 

24 

6.45E-04 

9.96E-04 

9.61E-04 

1.02E-03 

9.05E-04 

0.193195 

8.71E-04 

32 

9.09E-04 

7.73E-04 

6.86E-04 

5.98E-04 

7.42E-04 

0.178653 

8.35E-04 

40 

8.25E-04 

6.81E-04 

1.05E-03 

7.74E-04 

8.33E-04 

0.188349 

8.51E-04 

48 

9.20E-04 

8.44E-04 

9.98E-04 

7.25E-04 

8.72E-04 

0.133392 

8.05E-04 

56 

1.14E-03 

7.31E-04 

8.13E-04 

7.98E-04 

8.70E-04 

0.207804 

7.97E-04 

304 

4.09E-04 

4.38E-04 

7.67E-04 

6.22E-04 

5.59E-04 

0.091557 

5.28E-04 

312 

5.39E-04 

5.21E-04 

6.71E-04 

3.39E-04 

5.17E-04 

0.300373 

4.93E-04 

320 

4.43E-04 

5.03E-04 

4.66E-04 

3.90E-04 

4.51E-04 

0.263403 

4.86E-04 

328 

4.28E-04 

3.52E-04 

3.38E-04 

7.57E-04 

4.69E-04 

0.104947 

5.08E-04 

336 

3.25E-04 

2.65E-04 

5.42E-04 

4.32E-04 

3.91E-04 

0.418887 

4.90E-04 

344 

5.83E-04 

3.97E-04 

4.21E-04 

5.17E-04 

4.79E-04 

0.311981 

5.08E-04 

352 

4.09E-04 

4.38E-04 

7.67E-04 

6.22E-04 

5.59E-04 

0.179697 

5.28E-04 

Table  3b.  Details  of  PF  from  SS/SPT  method  calculated  at  8  sec  intervals 


Time 

(sec) 

SS/SPT 

MCS 

Run  1 

Run  2 

Run  3 

Run  4 

Mean 

C.O.V. 

8 

0.00652 

0.00764 

0.00616 

0.00529 

0.006402 

0.151619 

0.007024 

16 

0.01415 

0.01423 

0.01325 

0.01194 

0.013393 

0.079600 

0.013917 

24 

0.01924 

0.02208 

0.02084 

0.01998 

0.020535 

0.059432 

0.020785 

32 

0.02637 

0.02813 

0.02621 

0.02467 

0.026347 

0.053640 

0.027326 

40 

0.03280 

0.03343 

0.03439 

0.03071 

0.032834 

0.047451 

0.033949 

48 

0.03992 

0.03995 

0.04211 

0.03634 

0.03958 

0.060412 

0.040174 

56 

0.04864 

0.04557 

0.04834 

0.04249 

0.04626 

0.061983 

0.046292 

304 

0.18609 

0.17862 

0.18231 

0.17817 

0.181296 

0.020390 

0.183273 

312 

0.18875 

0.18150 

0.18733 

0.18226 

0.184959 

0.019571 

0.186726 

320 

0.19225 

0.18491 

0.19169 

0.18448 

0.188331 

0.022360 

0.189936 

328 

0.19511 

0.18819 

0.19470 

0.18702 

0.191256 

0.022200 

0.193087 

336 

0.19787 

0.19048 

0.19688 

0.19195 

0.194292 

0.018684 

0.196365 

344 

0.19995 

0.19219 

0.20036 

0.19474 

0.196812 

0.020341 

0.199516 

352 

0.20368 

0.19476 

0.20306 

0.19807 

0.199891 

0.021240 

0.202766 
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5.  SUMMARY,  CONCLUSIONS  AND  FUTURE 
WORK 

We  presented  an  improved  subset  simulation  with 
splitting  approach  to  calculate  the  time-dependent 
probability  of  failure  for  dynamic  systems  with  uncertain 
parameters  subjected  to  stochastic  excitation.  The  method 
partitions  the  high  dimensional  random  process  of  the 
response  into  a  series  of  correlated,  short  duration,  low 
dimensional  random  processes.  Subset  simulation  reduces 
the  computational  cost  by  introducing  appropriate 
intermediate  failure  sub-domains  to  express  the  low 
failure  probability  as  a  product  of  larger  conditional 
failure  probabilities.  Splitting  provides  an  efficient 
sampling  scheme  to  estimate  the  conditional  probabilities. 
The  proposed  subset  simulation  with  splitting  not  only 
estimates  the  time-dependent  probability  of  failure  at  a 
given  time  but  also  estimates  the  failure  rate  and 
cumulative  distribution  function  up  to  that  time  with 
approximately  the  same  cost.  We  used  the  vibratory 
response  of  a  vehicle  over  a  stochastic  terrain  to 
demonstrate  the  accuracy  and  efficiency  of  the  proposed 
approach. 

In  future  work,  we  plan  to  enhance  the  developed 
method  by  using  splitting  not  only  in  the  threshold 
direction  but  also  in  the  time  direction. 
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